rm(list = ls())
require(miceadds)
require(stats4)
require(lmtest)
require(bbmle)
require(msm)
require(nlWaldTest)
require(sandwich)
require(stargazer)
require(mfx)
require("multiwayvcov")
require(car)
require(alr3)
 
#data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
#data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
#data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)

data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)


dataframe2old=data.frame(data2old)
dataframe13=data.frame(data13)
dataframe15=data.frame(data15)


priors=c(0.5,0.6,0.75,0.85)

expframe=rbind(dataframe13,dataframe15,dataframe2old)
expframe=subset(expframe,uid>200)

#note that player 234 did not finish and so has no real data
expframe=subset(expframe,uid!=234)


#definitions
expframe$priorA=(expframe$qid==1)*0.5+(expframe$qid==5)*0.6+(expframe$qid==6)*0.75+(expframe$qid==7)*0.85
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==6)
expframe$Bresponse=(expframe$response==7)
expframe$stateA=(expframe$state==4)
expframe$stateB=(expframe$state==10)
expframe$correct=expframe$stateA*expframe$Aresponse+expframe$stateB*expframe$Bresponse
uidlist=unique(expframe$uid)

#find N
#**
length(uidlist)
#**

#fancy regression for easy clustering by Markus Conrad https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/
ols <- function(form, data, robust=FALSE, cluster=NULL,digits=3){
  r1 <- lm(form, data)
  if(length(cluster)!=0){
    data <- na.omit(data[,c(colnames(r1$model),cluster)])
    r1 <- lm(form, data)
  }
  X <- model.matrix(r1)
  n <- dim(X)[1]
  k <- dim(X)[2]
  if(robust==FALSE & length(cluster)==0){
    se <- sqrt(diag(solve(crossprod(X)) * as.numeric(crossprod(resid(r1))/(n-k))))
    res <- cbind(coef(r1),se)
  }
  if(robust==TRUE){
    u <- matrix(resid(r1))
    meat1 <- t(X) %*% diag(diag(crossprod(t(u)))) %*% X
    dfc <- n/(n-k)    
    se <- sqrt(dfc*diag(solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))))
    vcov <- solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))
    res <- cbind(coef(r1),se)
    }
  if(length(cluster)!=0){
    clus <- cbind(X,data[,cluster],resid(r1))
    colnames(clus)[(dim(clus)[2]-1):dim(clus)[2]] <- c(cluster,"resid")
    m <- dim(table(clus[,cluster]))
    dfc <- (m/(m-1))*((n-1)/(n-k))
    uclust  <- apply(resid(r1)*X,2, function(x) tapply(x, clus[,cluster], sum))
    se <- sqrt(diag(solve(crossprod(X)) %*% (t(uclust) %*% uclust) %*% solve(crossprod(X)))*dfc)   
    vcov <- solve(crossprod(X)) %*% (t(uclust) %*% uclust) %*% solve(crossprod(X))
    res <- cbind(coef(r1),se)
    }
  res <- cbind(res,res[,1]/res[,2],(1-pnorm(abs(res[,1]/res[,2])))*2)
  res1 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res)),nrow=dim(res)[1])
  rownames(res1) <- rownames(res)
  colnames(res1) <- c("Estimate","Std. Error","t value","Pr(>|t|)")
  out=list(res1,vcov)
  return(out)
}




#Find thresholds and look for dropping

priorchangeframe=data.frame(uid=vector('numeric'),pA0.5=vector('numeric'),pA0.5lower=vector('numeric'),pA0.5upper=vector('numeric'),pA0.6=vector('numeric'),pA0.75=vector('numeric'),pA0.85=vector('numeric'),ChangePost0.5_0.6=vector('numeric'),Chisq0.5_0.6=vector('numeric'),Prob0.5_0.6=vector('numeric'),ChangePost0.5_0.75=vector('numeric'),Chisq0.5_0.75=vector('numeric'),Prob0.5_0.75=vector('numeric'),ChangePost0.5_0.85=vector('numeric'),Chisq0.5_0.85=vector('numeric'),Prob0.5_0.85=vector('numeric'))

for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
model=lm(Aresponse~factor(stateA+priorA)+0,data=playframe)
na1=model$coef[1]
na2=model$coef[2]
na3=model$coef[3]
na4=model$coef[4]
a1=model$coef[5]
a2=model$coef[6]
a3=model$coef[7]
a4=model$coef[8]

if(sum(playframe$Aresponse*(playframe$priorA==0.6))==0){
postAhigh0.6=0

}else{
postAhigh0.6=(0.6*a2)/(0.6*a2+0.4*na2)
}


if(sum(playframe$Aresponse*(playframe$priorA==0.75))==0){
postAhigh0.75=0
}else{
postAhigh0.75=(0.75*a3)/(0.75*a3+0.25*na3)
}

if(sum(playframe$Aresponse*(playframe$priorA==0.85))==0){
postAhigh0.85=0
}else{
postAhigh0.85=(0.85*a4)/(0.85*a4+0.15*na4)
}

if(sum(playframe$Aresponse*(playframe$priorA==0.5))==0){

postAhigh0.5=0
postAhigh0.5upper=0
postAhigh0.5lower=0

diff0.5_0.6=postAhigh0.5-postAhigh0.6
chisq0.5_0.6=1e7*(postAhigh0.5-postAhigh0.6>0)
diff0.5_0.6p=7.888609e-31*(postAhigh0.5-postAhigh0.6>0)+(1-7.888609e-31)*(postAhigh0.5-postAhigh0.6==0)

diff0.5_0.75=postAhigh0.5-postAhigh0.75
chisq0.5_0.75=1e7*(postAhigh0.5-postAhigh0.75>0)
diff0.5_0.75p=7.888609e-31*(postAhigh0.5-postAhigh0.75>0)+(1-7.888609e-31)*(postAhigh0.5-postAhigh0.75==0)

diff0.5_0.85=postAhigh0.5-postAhigh0.85
chisq0.5_0.85=1e7*(postAhigh0.5-postAhigh0.85>0)
diff0.5_0.85p=7.888609e-31*(postAhigh0.5-postAhigh0.85>0)+(1-7.888609e-31)*(postAhigh0.5-postAhigh0.85==0)


}else{
postAhigh0.5=(0.5*a1)/(0.5*a1+0.5*na1)

vcovmod=vcovHC(model, type="HC3")

threshsterr=as.numeric(deltamethod(~(0.5*x1)/(0.5*x1+0.5*x2),coef(model)[c(5,1)],vcovmod[c(5,1),c(5,1)]))
postAhigh0.5upper=postAhigh0.5+1.96*threshsterr
postAhigh0.5lower=postAhigh0.5-1.96*threshsterr

if(sum(playframe$Aresponse)==length(playframe$Aresponse)|sum(playframe$Aresponse)==0){
diff0.5_0.6=0
chisq0.5_0.6=0
diff0.5_0.6p=1

diff0.5_0.75=0
chisq0.5_0.75=0
diff0.5_0.75p=1

diff0.5_0.85=0
chisq0.5_0.85=0
diff0.5_0.85p=1

}else{

diff0.5_0.6=postAhigh0.5-postAhigh0.6
if (sum(playframe$Aresponse*(playframe$priorA==0.5))==sum(playframe$Aresponse*(playframe$priorA==0.6))){
chisq0.5_0.6=0
diff0.5_0.6p=1
}else{
test0.5_0.6=nlWaldtest(model,"(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.6*b[6])/(0.6*b[6]+0.4*b[2])",Vcov=vcovmod)
chisq0.5_0.6=as.numeric(test0.5_0.6[1])
diff0.5_0.6p=as.numeric(test0.5_0.6[4])
}

diff0.5_0.75=postAhigh0.5-postAhigh0.75
test0.5_0.75=nlWaldtest(model,"(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.75*b[7])/(0.75*b[7]+0.25*b[3])",Vcov=vcovmod)
chisq0.5_0.75=as.numeric(test0.5_0.75[1])
diff0.5_0.75p=as.numeric(test0.5_0.75[4])

diff0.5_0.85=postAhigh0.5-postAhigh0.85
test0.5_0.85=nlWaldtest(model,"(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.85*b[8])/(0.85*b[8]+0.15*b[4])",Vcov=vcovmod)
chisq0.5_0.85=as.numeric(test0.5_0.85[1])
diff0.5_0.85p=as.numeric(test0.5_0.85[4])

}

}
addframe=data.frame(uid=uidlist[i],pA0.5=postAhigh0.5,pA0.5lower=postAhigh0.5lower,pA0.5upper=postAhigh0.5upper,pA0.6=postAhigh0.6,pA0.75=postAhigh0.75,pA0.85=postAhigh0.85,ChangePost0.5_0.6=diff0.5_0.6,Chisq0.5_0.6=chisq0.5_0.6,Prob0.5_0.6=diff0.5_0.6p,ChangePost0.5_0.75=diff0.5_0.75,Chisq0.5_0.75=chisq0.5_0.75,Prob0.5_0.75=diff0.5_0.75p,ChangePost0.5_0.85=diff0.5_0.85,Chisq0.5_0.85=chisq0.5_0.85,Prob0.5_0.85=diff0.5_0.85p)


priorchangeframe=rbind(priorchangeframe,addframe)
}

percent=vector('numeric')
percent[1]=100*sum((priorchangeframe$pA0.5<0.6))/length(priorchangeframe$uid)
percent[2]=100*sum((priorchangeframe$pA0.5>=0.6)*(priorchangeframe$pA0.5<0.75))/length(priorchangeframe$uid)
percent[3]=100*sum((priorchangeframe$pA0.5>=0.75)*(priorchangeframe$pA0.5<0.85))/length(priorchangeframe$uid)
percent[4]=100*sum((priorchangeframe$pA0.5>=0.85)*(priorchangeframe$pA0.5<=1))/length(priorchangeframe$uid)
thresholds=c("[0.5,0.6)","[0.6,0.75)","[0.75,0.85)","[0.85,1]")

thresholdsframe=data.frame(Theshold=thresholds,Percent=percent)

#***********************************
stargazer(thresholdsframe,rownames=FALSE,summary=FALSE)
#***********************************


#Make some graphs of posteriors with subjects that have different thresholds
above0.6=subset(priorchangeframe, pA0.5lower>0.6)

uidlistabove0.6=unique(above0.6$uid)
above0.6frame0.5=subset(expframe,priorA==0.5&uid %in%uidlistabove0.6)
above0.6frame0.6=subset(expframe,priorA==0.6&uid %in%uidlistabove0.6)
#for p values

comboframe0.6=rbind(above0.6frame0.5,above0.6frame0.6)
comboframe0.6$highpri=(comboframe0.6$priorA==0.6)
diffmodelA=ols(stateB~Aresponse+highpri+highpri*Aresponse,data=comboframe0.6,robust=TRUE, cluster="uid")
diffmodelB=ols(stateA~Bresponse+highpri+highpri*Bresponse,data=comboframe0.6,robust=TRUE, cluster="uid")

#number of significant violations
priorchangeframeabove0.6=subset(priorchangeframe,uid %in% uidlistabove0.6)

#****************************************
sum(priorchangeframeabove0.6$Prob0.5_0.6<0.05)/length(priorchangeframeabove0.6$Prob0.5_0.6)
#****************************************

#for graphs
above0.6frame=subset(expframe,uid %in% uidlistabove0.6)
model0.6=lm(Aresponse~factor(stateA+priorA)+0, data=above0.6frame)
vcovmod0.6cluster=cluster.vcov(model0.6, above0.6frame$uid)

correctfrac0.5A=0.5*model0.6$coefficient[5]/(0.5*model0.6$coefficient[5]+0.5*model0.6$coefficient[1])
correctfrac0.5B=0.5*(1-model0.6$coefficient[1])/(0.5*(1-model0.6$coefficient[1])+0.5*(1-model0.6$coefficient[5]))
correctfrac0.6A=0.6*model0.6$coefficient[6]/(0.6*model0.6$coefficient[6]+0.4*model0.6$coefficient[2])
correctfrac0.6B=0.4*(1-model0.6$coefficient[2])/(0.4*(1-model0.6$coefficient[2])+0.6*(1-model0.6$coefficient[6]))
correctfrac0.5=c(correctfrac0.5A,correctfrac0.5B)
correctfrac0.6=c(correctfrac0.6A,correctfrac0.6B)

sterr=vector("numeric")

sterr[1]=as.numeric(deltamethod(~0.5*x5/(0.5*x5+0.5*x1),model0.6$coefficients,vcovmod0.6cluster))
sterr[2]=as.numeric(deltamethod(~0.5*(1-x1)/(0.5*(1-x1)+0.5*(1-x5)),model0.6$coefficients,vcovmod0.6cluster))
sterr[3]=as.numeric(deltamethod(~0.6*x6/(0.6*x6+0.4*x2),model0.6$coefficients,vcovmod0.6cluster))
sterr[4]=as.numeric(deltamethod(~0.4*(1-x2)/(0.4*(1-x2)+0.6*(1-x6)),model0.6$coefficients,vcovmod0.6cluster))

correctfrac0.5high=correctfrac0.5+c(sterr[1],sterr[2])
correctfrac0.6high=correctfrac0.6+c(sterr[3],sterr[4])
correctfrac0.5low=correctfrac0.5-c(sterr[1],sterr[2])
correctfrac0.6low=correctfrac0.6-c(sterr[3],sterr[4])

#*****************************
barCenters=barplot(matrix(c(rbind(correctfrac0.5,correctfrac0.6)),nr=2), beside=T, 
        col=c("blue","orange"), 
        names.arg=c("P(1|A)","P(2|B)"),ylim = c(0, 1),xpd = FALSE,cex.axis=2.5,cex.names=2.5)
 par(xpd=TRUE)
legend(0.7,-0.06,ncol=3, c("Prior","0.5","0.6"), pch=15, 
       col=c("white","blue","orange"), 
       bty="n",cex=2.5)

box(bty="l")

segments(barCenters,c(rbind(correctfrac0.5low,correctfrac0.6low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.6high)), lwd = 1.5)

arrows(barCenters, c(rbind(correctfrac0.5low,correctfrac0.6low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.6high)), lwd = 1.5, angle = 90,
       code = 3, length = 0.05)

#*******************************

#Difference p-values
Textstr=vector("character")
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.6*b[6])/(0.6*b[6]+0.4*b[2])"
nlWaldtest(model0.6,texts=Textstr,Vcov=vcovmod0.6cluster)

Textstr[1]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.6*(1-b[2]))/(0.6*(1-b[2])+0.4*(1-b[6]))"
nlWaldtest(model0.6,texts=Textstr,Vcov=vcovmod0.6cluster)

#Next Graph
above0.75=subset(priorchangeframe, pA0.5lower>0.75)
uidlistabove0.75=unique(above0.75$uid)
above0.75frame0.5=subset(expframe,priorA==0.5&uid %in%uidlistabove0.75)
above0.75frame0.75=subset(expframe,priorA==0.75&uid %in%uidlistabove0.75)
#for p values

comboframe0.75=rbind(above0.75frame0.5,above0.75frame0.75)
comboframe0.75$highpri=(comboframe0.75$priorA==0.75)
diffmodelA=ols(stateB~Aresponse+highpri+highpri*Aresponse,data=comboframe0.75,robust=TRUE, cluster="uid")
diffmodelB=ols(stateA~Bresponse+highpri+highpri*Bresponse,data=comboframe0.75,robust=TRUE, cluster="uid")

#number of significant violations
priorchangeframeabove0.75=subset(priorchangeframe,uid %in% uidlistabove0.75)
#****************************************
sum(priorchangeframeabove0.75$Prob0.5_0.75<0.05)/length(priorchangeframeabove0.75$Prob0.5_0.75)
#****************************************


above0.75=subset(priorchangeframe, pA0.5lower>0.75)
uidlistabove0.75=unique(above0.75$uid)


#for graphs
above0.75frame=subset(expframe,uid %in% uidlistabove0.75)
model0.75=lm(Aresponse~factor(stateA+priorA)+0, data=above0.75frame)
vcovmod0.75cluster=cluster.vcov(model0.75, above0.75frame$uid)

correctfrac0.5A=0.5*model0.75$coefficient[5]/(0.5*model0.75$coefficient[5]+0.5*model0.75$coefficient[1])
correctfrac0.5B=0.5*(1-model0.75$coefficient[1])/(0.5*(1-model0.75$coefficient[1])+0.5*(1-model0.75$coefficient[5]))
correctfrac0.75A=0.75*model0.75$coefficient[7]/(0.75*model0.75$coefficient[7]+0.25*model0.75$coefficient[3])
correctfrac0.75B=0.25*(1-model0.75$coefficient[3])/(0.25*(1-model0.75$coefficient[3])+0.75*(1-model0.75$coefficient[7]))
correctfrac0.5=c(correctfrac0.5A,correctfrac0.5B)
correctfrac0.75=c(correctfrac0.75A,correctfrac0.75B)

sterr=vector("numeric")

sterr[1]=as.numeric(deltamethod(~0.5*x5/(0.5*x5+0.5*x1),model0.75$coefficients,vcovmod0.75cluster))
sterr[2]=as.numeric(deltamethod(~0.5*(1-x1)/(0.5*(1-x1)+0.5*(1-x5)),model0.75$coefficients,vcovmod0.75cluster))
sterr[3]=as.numeric(deltamethod(~0.75*x7/(0.75*x7+0.25*x3),model0.75$coefficients,vcovmod0.75cluster))
sterr[4]=as.numeric(deltamethod(~0.25*(1-x3)/(0.25*(1-x3)+0.75*(1-x7)),model0.75$coefficients,vcovmod0.75cluster))

correctfrac0.5high=correctfrac0.5+c(sterr[1],sterr[2])
correctfrac0.75high=correctfrac0.75+c(sterr[3],sterr[4])
correctfrac0.5low=correctfrac0.5-c(sterr[1],sterr[2])
correctfrac0.75low=correctfrac0.75-c(sterr[3],sterr[4])

#*****************************
barCenters=barplot(matrix(c(rbind(correctfrac0.5,correctfrac0.75)),nr=2), beside=T, 
        col=c("blue","orange"), 
        names.arg=c("P(1|A)","P(2|B)"),ylim = c(0, 1),xpd = FALSE,cex.axis=2.5,cex.names=2.5)
 par(xpd=TRUE)
legend(0.7,-0.06,ncol=3, c("Prior","0.5","0.75"), pch=15, 
       col=c("white","blue","orange"), 
       bty="n",cex=2.5)
box(bty="l")

segments(barCenters,c(rbind(correctfrac0.5low,correctfrac0.75low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.75high)), lwd = 1.5)

arrows(barCenters, c(rbind(correctfrac0.5low,correctfrac0.75low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.75high)), lwd = 1.5, angle = 90,
       code = 3, length = 0.05)

#*******************************
 
#Difference p-values
Textstr=vector("character")
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.75*b[7])/(0.75*b[7]+0.25*b[3])"
nlWaldtest(model0.75,texts=Textstr,Vcov=vcovmod0.75cluster)

Textstr[1]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.75*(1-b[3]))/(0.75*(1-b[7])+0.25*(1-b[7]))"
nlWaldtest(model0.75,texts=Textstr,Vcov=vcovmod0.75cluster)


#Next Graph
above0.85=subset(priorchangeframe, pA0.5lower>0.85)
uidlistabove0.85=unique(above0.85$uid)
above0.85frame0.5=subset(expframe,priorA==0.5&uid %in%uidlistabove0.85)
above0.85frame0.85=subset(expframe,priorA==0.85&uid %in%uidlistabove0.85)

#for p values

comboframe0.85=rbind(above0.85frame0.5,above0.85frame0.85)
comboframe0.85$highpri=(comboframe0.85$priorA==0.85)
diffmodelA=ols(stateB~Aresponse+highpri+highpri*Aresponse,data=comboframe0.85,robust=TRUE, cluster="uid")
diffmodelB=ols(stateA~Bresponse+highpri+highpri*Bresponse,data=comboframe0.85,robust=TRUE, cluster="uid")

#number of significant violations
priorchangeframeabove0.85=subset(priorchangeframe,uid %in% uidlistabove0.85)

#****************************************
sum(priorchangeframeabove0.85$Prob0.5_0.85<0.05)/length(priorchangeframeabove0.85$Prob0.5_0.85)
#****************************************

#for graphs

above0.85=subset(priorchangeframe, pA0.5lower>0.85)
uidlistabove0.85=unique(above0.85$uid)


#for graphs
above0.85frame=subset(expframe,uid %in% uidlistabove0.85)
model0.85=lm(Aresponse~factor(stateA+priorA)+0, data=above0.85frame)
vcovmod0.85cluster=cluster.vcov(model0.85, above0.85frame$uid)

correctfrac0.5A=0.5*model0.85$coefficient[5]/(0.5*model0.85$coefficient[5]+0.5*model0.85$coefficient[1])
correctfrac0.5B=0.5*(1-model0.85$coefficient[1])/(0.5*(1-model0.85$coefficient[1])+0.5*(1-model0.85$coefficient[5]))
correctfrac0.85A=0.85*model0.85$coefficient[8]/(0.85*model0.85$coefficient[8]+0.15*model0.85$coefficient[4])
correctfrac0.85B=0.15*(1-model0.85$coefficient[4])/(0.15*(1-model0.85$coefficient[4])+0.85*(1-model0.85$coefficient[8]))
correctfrac0.5=c(correctfrac0.5A,correctfrac0.5B)
correctfrac0.85=c(correctfrac0.85A,correctfrac0.85B)

sterr=vector("numeric")

sterr[1]=as.numeric(deltamethod(~0.5*x5/(0.5*x5+0.5*x1),model0.85$coefficients,vcovmod0.85cluster))
sterr[2]=as.numeric(deltamethod(~0.5*(1-x1)/(0.5*(1-x1)+0.5*(1-x5)),model0.85$coefficients,vcovmod0.85cluster))
sterr[3]=as.numeric(deltamethod(~0.85*x8/(0.85*x8+0.15*x4),model0.85$coefficients,vcovmod0.85cluster))
sterr[4]=as.numeric(deltamethod(~0.15*(1-x4)/(0.15*(1-x4)+0.85*(1-x8)),model0.85$coefficients,vcovmod0.85cluster))

correctfrac0.5high=correctfrac0.5+c(sterr[1],sterr[2])
correctfrac0.85high=correctfrac0.85+c(sterr[3],sterr[4])
correctfrac0.5low=correctfrac0.5-c(sterr[1],sterr[2])
correctfrac0.85low=correctfrac0.85-c(sterr[3],sterr[4])

#*****************************
barCenters=barplot(matrix(c(rbind(correctfrac0.5,correctfrac0.85)),nr=2), beside=T, 
        col=c("blue","orange"), 
        names.arg=c("P(1|A)","P(2|B)"),ylim = c(0, 1),xpd = FALSE,cex.axis=2.5,cex.names=2.5)
 par(xpd=TRUE)
legend(0.7,-0.06,ncol=3, c("Prior","0.5","0.85"), pch=15, 
       col=c("white","blue","orange"), 
       bty="n",cex=2.5)

box(bty="l")

segments(barCenters,c(rbind(correctfrac0.5low,correctfrac0.85low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.85high)), lwd = 1.5)

arrows(barCenters, c(rbind(correctfrac0.5low,correctfrac0.85low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.85high)), lwd = 1.5, angle = 90,
       code = 3, length = 0.05)

#*******************************

#Difference p-values
Textstr=vector("character")
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.85*b[8])/(0.85*b[8]+0.15*b[4])"
nlWaldtest(model0.85,texts=Textstr,Vcov=vcovmod0.85cluster)

Textstr[1]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.85*(1-b[4]))/(0.85*(1-b[4])+0.15*(1-b[8]))"
nlWaldtest(model0.85,texts=Textstr,Vcov=vcovmod0.85cluster)


############################Aggregate LIP Test#######################################
above0.85frameall=subset(expframe,uid %in% uidlistabove0.85 )
modelLIP=lm(Aresponse~factor(stateA+priorA)+0,data=above0.85frameall)


vcovmodLIP=vcovHC(modelLIP, type="HC3")
vcovmodLIPcluster=cluster.vcov(modelLIP, above0.85frameall$uid)

Textstr=vector("character")
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.6*b[6])/(0.6*b[6]+0.4*b[2])"
Textstr[2]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.6*(1-b[2]))/(0.6*(1-b[2])+0.4*(1-b[6]))"
Textstr[3]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.75*b[7])/(0.75*b[7]+0.25*b[3])"
Textstr[4]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.75*(1-b[3]))/(0.75*(1-b[3])+0.25*(1-b[7]))"
Textstr[5]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.85*b[8])/(0.85*b[8]+0.15*b[4])"
Textstr[6]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.85*(1-b[4]))/(0.85*(1-b[4])+0.15*(1-b[8]))"

#nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIP)
waldmodel=nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)
waldmodel
nlWaldtest(modelLIP,texts=Textstr)




Textstr=vector("character")
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.6*b[6])/(0.6*b[6]+0.4*b[2])"
nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)
Textstr[1]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.6*(1-b[2]))/(0.6*(1-b[2])+0.4*(1-b[6]))"
nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.75*b[7])/(0.75*b[7]+0.25*b[3])"
nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)
Textstr[1]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.75*(1-b[3]))/(0.75*(1-b[3])+0.25*(1-b[7]))"
nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)
Textstr[1]="(0.5*b[5])/(0.5*b[5]+0.5*b[1])-(0.85*b[8])/(0.85*b[8]+0.15*b[4])"
nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)
Textstr[1]="(0.5*(1-b[1]))/(0.5*(1-b[1])+0.5*(1-b[5]))-(0.85*(1-b[4]))/(0.85*(1-b[4])+0.15*(1-b[8]))"
nlWaldtest(modelLIP,texts=Textstr,Vcov=vcovmodLIPcluster)




